#include <iostream>
#include <cmath>
#include <Eigen/Dense>
#include <ctime>
#include <vector>
#include <numeric>
#include <fstream>

using namespace std;
using namespace Eigen;


template <typename T>
struct ColComparator
{
	const MatrixBase<T>& matrix;
	const int row;

	ColComparator(MatrixBase<T>& m, int row)
		: matrix(m.derived()), row(row)
	{
	}

	bool operator()(int i, int j) const
	{
		return matrix(row, i) < matrix(row, j);
	}
};
//need verify
void sortByCol(MatrixXd& m, int n)
{
	MatrixXd tmp = m;
	VectorXi indice = VectorXi::LinSpaced(m.cols(), 0, m.cols() - 1);
	std::sort(indice.data(), indice.data() + indice.size(),
		ColComparator<MatrixXd>(m, n));
	for (int i = 0; i < indice.size(); i++)
		m.col(i) = tmp.col(indice(i));
}

//need verify
void premuteByCol(MatrixXd& m)
{
	MatrixXd tmp = m;
	VectorXi indice = VectorXi::LinSpaced(m.cols(), 0, m.cols() - 1);
	std::random_shuffle(indice.data(), indice.data() + indice.size());
	for (int i = 0; i < indice.size(); i++)
		m.col(i) = tmp.col(indice(i));
}

MatrixXd* readFromFile(const string& fileName)
{
	ifstream ifd(fileName, ifstream::in);
	string line;
	stringstream ss;
	int rows, cols;
	if (!ifd.is_open())
	{
		cerr << "Open Error: " << fileName;
		exit(1);
	}
	getline(ifd, line);
	ss << line;
	ss >> rows >> cols;
	ss.clear();
	MatrixXd* m = new MatrixXd(rows, cols);
	int c = 0;
	while (std::getline(ifd, line))
	{
		ss << line;
		for (int i = 0; i < cols; ++i)
		{
			ss >> (*m)(c, i);
		}
		ss.clear();
		++c;
	}
	ifd.close();
	m->transposeInPlace();
	return m;
}

struct Sigmoid
{
	double operator()(double x) const
	{
		return 1 / (1 + exp(-x));
	}
};
//nh: number of elements in Hide Layer
//xi: for stop iterating
//eta: for train ratio
VectorXd train(int nh, const VectorXd &X, const VectorXd &Y, MatrixXd &V, MatrixXd &W, VectorXd &_theta, VectorXd &_gamma, double eta)
{
	VectorXd B(nh);
	VectorXd Y_bar(Y.size());
	B = V*X - _gamma;
	B = B.unaryExpr(Sigmoid()).eval();//hidden layer 
	Y_bar = W*B - _theta;
	Y_bar = Y_bar.unaryExpr(Sigmoid()).eval();//output layer

	VectorXd G = Y_bar.array()*(1 - Y_bar.array())*(Y.array() - Y_bar.array());
	VectorXd E = B.array()*(1 - B.array())*(W.transpose()*G).array();

	W = (W + eta*G*B.transpose()).eval();
	_theta = (_theta - eta*G).eval();
	V = (V + eta*E*X.transpose()).eval();
	_gamma = (_gamma - eta*E);
	Y_bar = W*B - _theta;
	Y_bar = Y_bar.unaryExpr(Sigmoid()).eval();
	return Y_bar;
}

void bp(int nh, const MatrixXd &train_data, double xi, double eta, MatrixXd &V, MatrixXd &W, VectorXd &_theta, VectorXd &_gamma)
{

	//prepare data
	MatrixXd X = train_data.block(0, 0, train_data.rows() - 1, train_data.cols());
	VectorXi labels = train_data.block(train_data.rows() - 1, 0, 1, train_data.cols()).transpose().cast<int>();
	int maxElem = labels.maxCoeff();
	MatrixXd Y = MatrixXd::Zero(maxElem, labels.size());

	for (int i = 0; i < labels.size(); ++i)
		Y(labels(i) - 1, i) = 1;
	vector<double> eList;

	V = MatrixXd::Random(nh, X.rows()); //weight for input layer
	W = MatrixXd::Random(Y.rows(), nh); //weight for hide layer
	_theta = VectorXd::Random(Y.rows()); //threshold for hide layer
	_gamma = VectorXd::Random(nh);//threshold for output layer

	double avg_ek;//average of Ek
	double pre = 0, tmp;//for stop
	do
	{
		eList.clear();
		for (auto i = 0; i < X.cols(); ++i)
		{
			VectorXd &Y_bar = train(nh, X.col(i), Y.col(i), V, W, _theta, _gamma, eta);
			eList.push_back((Y.col(i) - Y_bar).squaredNorm());
		}
		avg_ek = 0;
		for (auto p = eList.cbegin(); p != eList.cend(); ++p)
			avg_ek += *p;
		avg_ek /= eList.size();
		tmp = abs(avg_ek - pre);
		pre = avg_ek;
	} while (tmp > xi);
}

double validate(const MatrixXd &test_data, int clzNum, const MatrixXd &V, const MatrixXd &W, const VectorXd& _theta, const VectorXd &_gamma)
{
	MatrixXd T = test_data.block(0, 0, test_data.rows() - 1, test_data.cols());
	VectorXi labels = test_data.block(test_data.rows() - 1, 0, 1, test_data.cols()).transpose().cast<int>();
	int maxElem = labels.maxCoeff();
	MatrixXd Labels = MatrixXd::Identity(clzNum, clzNum);
	int count = 0;
	for (int i = 0; i < T.cols(); ++i)
	{
		VectorXd B(V.rows());
		VectorXd Y_bar;
		B = V*T.col(i) - _gamma;
		B = B.unaryExpr(Sigmoid()).eval();//hidden layer 
		Y_bar = W*B - _theta;
		Y_bar = Y_bar.unaryExpr(Sigmoid()).eval();//output layer

		VectorXd err(clzNum);
		for (int j = 0; j < clzNum; ++j)
			err(j) = (Y_bar - Labels.col(j)).squaredNorm();
		int maxIndex;
		err.minCoeff(&maxIndex);
		if (maxIndex == (labels(i) - 1))
			count++;
	}
	return (double)count / test_data.cols();
}

int main()
{
	MatrixXd* d;
	std::srand(unsigned(std::time(0)));
	d = readFromFile("iris.data");
	double rec = 0;
	for (int i = 0; i < 100; i++)
	{
		premuteByCol(*d);
		MatrixXd train_data, verify, test;
		train_data = d->block(0, 0, d->rows(), 20);
		test = d->block(0, 20, d->rows(), 130);

		MatrixXd V, W;
		VectorXd _theta, _gamma;
		bp(4, train_data, 0.00001, 0.1, V, W, _theta, _gamma);
		int clzNum = d->row(d->rows() - 1).maxCoeff();
		double tmp = validate(test, clzNum, V, W, _theta, _gamma);
		cout << tmp << endl;
		rec += tmp;
	}
	cout << rec / 100 << endl;
	return 0;
}
